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Abstract 

We study the Gierer-Meinhardt model of reaction-diffusion on a 
site-disordered square lattice. Let p be the site occupation probability 
of the square lattice. For p greater than a critical value p c , the steady 
state consists of stripe-like patterns with long-range connectivity. For 
p < Pc, the connectivity is lost. The value of p c is found to be 
much greater than that of the site percolation threshold for the square 
lattice. In the vicinity of p c , the cluster-related quantities exhibit 
power-law scaling behaviour. The method of finite-size scaling is used 
to determine the values of the fractal dimension df, the ratio, ^, of the 
average cluster size exponent 7 and the correlation length exponent u 
and also v itself. The values appear to indicate that the disordered 
GM model belongs to the universality class of ordinary percolation. 

PACS number(s): 05.70. Ln 

I. Introduction 

Phase transitions in the non-equilibrium are not as well understood as in 
the case of equilibrium systems. The percolation model provides an exam- 
ple of phase transition in disordered systems pi The model is simply 
defined. Consider a square lattice of sites. The lattice can be made disor- 
dered by assuming that each site is present with probability p and absent 
with probability 1 — p. There is a critical value p c of p, called the perco- 
lation threshold, below which long-range connectivity in the system under 
study is missing and above which the connectivity is present with probability 
one. The percolation transition is analogous to equilibrium thermodynamic 
phase transitions and exhibits "critical" phenomena in the vicinity of the 
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percolation threshold. Percolation-like transitions have been observed in the 
non-equilibrium in a surface-reaction model ||, an autocatalytic reaction 
model H and chemical wave propagation on a lattice of excitable and non- 
excitable clusters ||. Ziff, Gulari and Barshad(ZGB) || have studied a 
kinetic reaction model describing the chemical reaction CO+O — >C0 2 on a 
catalytic surface. The surface is represented by a square lattice, the sites of 
which can be singly occupied by a CO molecule or an O atom, let X C o an d 
1-Xco be the molecular concentrations of CO and O2 molecules in the gas 
above the surface. CO and O2 molecules are added randomly to the lattice 
in the relative ratio of their molecular concentrations. Two adjacent empty 
sites are necessary for the adsorption of an oxygen molecule. The reaction 
consists of CO — O pairs being randomly selected from nearest-neighbour 
(n.n) lattice sites occupied by O atoms and CO molecules. The reaction 
produces C0 2 molecules which escape from the surface. ZGB found that 
when X C o is l ess than a certain value X%, the lattice is completely covered 
by oxygen atoms in the steady state and no further reaction takes place. For 
Xco greater than X2, the lattice is filled with CO. The system exists in a 
catalytically active mixed-phase only in the range X\ < Xco < X 2 - The 
transitions at Xi(X 2 ) have been found to be of second (first) order. Later 
studies II have shown that the second-order phase transition belongs to 
the universality class of Reggeon-field theory (RFT) to which directed perco- 
lation (DP) also belongs. Aukrust et al ||] have studied an irreversible kinetic 
reaction model for a one-component autocatalytic reaction A + A— > A 2 . In 
this model, an atom adsorbing on a lattice site reacts, with probablity 1-p, 
with one of the n.n. atoms, if present. After the reaction, the two atoms 
leave the lattice. Otherwise the atom occupies the site. As p is varied, the 
model undergoes a second-order kinetic phase transition from a partial occu- 
pation (chemically-reactive state) of the lattice to a completely covered inert 
state. Again, the phase transition has been shown to belongs to the univer- 
sality class of Reggeon-field theory-directed percolation. Sendina-Nadal et al 
have studied the propagation of chemical waves in a disordered excitable 
medium in terms of percolation theory. They performed an experiment with 
Belousov-Zhabotinsky reagent catalyzed by Ru{bpy) complex which is sensi- 
tive to visible light thus making it possible to control the excitability of the 
system experimentally. The system consists of excitable (black) and nonex- 
citable (white) clusters. The effective wave front velocity has been observed 
to jump from finite values at a threshold p = p c , where p is the proportion 
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of black sites. At p c , a cluster of black sites spans the medium. The data ob- 
tained appear to indicate that the percolation process is similar to ordinary 
percolation. In this paper, we study the Gierer-Meinhardt (GM) model || of 
reaction-diffusion (RD) in two dimensions (2d) and show that a percolation- 
like transition occurs in the non-equilibrium steady state as the underlying 
lattice is made disordered. In the vicinity of the percolation threshold p c , 
evidence of power-law scaling behaviour of cluster-related quantities is ob- 
tained. The percolation threshold appears to be in the universality class of 
ordinary percolation. The study is based on the finite-size scaling analysis of 
the percolation problem. In Section II, we define the GM model of RD and 
describe the various results obtained by us. Section III contains concluding 
remarks. 



II. Reaction-diffusion in the presence of disor- 
der 

Turing || made a remarkable suggestion, several years back, that diffusion 
need not always act to smooth out concentration differences in a chemical 
system. Two interacting chemicals can generate a stable, inhomogeneous 
pattern if one of the substances (the inhibitor) diffuses much faster than the 
other (the activator). The activator is autocatalytic, i.e, a small increase in 
its concentration 'a' over its homogeneous steady-state concentration leads to 
a further increse of 'a'. The activator besides promoting its own production 
also promotes the production of the inhibitor. The inhibitor, as the name 
implies, is antagonistic to the activator and inhibits its production. Suppose, 
originally the system is in a homogeneous steady state. A local increase in 
the activator concentration leads to a further increase in the concentration of 
the activator due to autocatalysis. The concentration of the inhibitor is also 
increased locally. The inhibitor, having a diffusion coefficient much larger 
than that of the activator, diffuses faster to the surrounding region and pre- 
vents the activator from coming there. This process of autocatalysis and 
long-rang inhibition finally lead to a stationary state consisting of islands of 
high activator concentration within a region of high inhibitor concentration. 
The islands constitute what is known as the Turing pattern. In some Turing 
systems, there is a possibility of saturation of the autocatalytic process. In 
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this case, the inhibitor production is limited. Locally activated regions now 
have activated neighbours. Nonactivated areas are, however, close by into 
which the inhibitor can diffuse. This results in the emergence of stripe-like 
patterns which have long-range connectivity. We study the effect of disorder 
on the connectivity. The Turing patterns are regions of chemical concentra- 
tion gradients. Turing's original idea was that the gradients are responsible 
for the patterns seen in biological systems. Later, several physical, chemi- 



cal and biological systems have been identified [|T0j, [11], |T|, |H| [14], [15|] which 
exhibit Turing patterns. 

The GM model of RD is described by the following two partial differential 
equations : 

da a 2 . 

m = DaAa + pa aTk^jh-^ a (la) 

dh 

— = D h Ah + p h a 2 - p h h (lb) 
at 

where A is the Laplacian given by A = + Jp- , 'a' and 'h' denote the 
concentrations of the activator and the inhibitor, D a , Dh, are the respec- 
tive diffusion coefficients, p a , ph are the removal rates and p a , Ph are the 
cross- react ion coefficients. The conditions for the formation of stable Turing 
patterns are D h » D a and ph > [§■ We also assume that j^- = p a 
and ph = Ph- In this case, the steady state solution of equations (la) and 
(lb) is given by (a,h) = (l,l), i.e., the steady state is homogeneous. We use 
the parameter values D a = 0.005, D h = 0.2, p h = Ph = 0.02, p a = 0.0125 
and k a = 0.25 for our study. Bose and Chaudhuri have studied the 
effect of disorder on the formation of Turing islands in the GM model with 
k a = 0. We extend this study to the case of k a ^ 0. As in Ref. [L(J, a very 



simple discretization scheme is used. The lattice chosen is of size LxL, with 
L ranging from 20 to 70 in steps of 10. Also, periodic boundary conditions 
are used. Disorder is introduced into the underlying substrate (lattice) by 
stipulating that the probability of a paticular site being part of the RD net- 
work is p. The disordered lattice configuration is generated with the help of 



a random number generator [T^|. For the disordered lattice, the Laplacian is 
written as 



A a(xij, t) = iocc(i + 1, j) x ( a(x i+ ij, t) - a(xij, t) ) 
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+ iocc(i, j + 1) x ( a(x ij+ i, t) - a(xij, t) ) 
+ iocc{i — 1, j) x ( a(xi-ij, t) — a(xij, t) ) 
+ iocc(i, j -1) x (a(%_i, £) - a(a^, t) ) (2) 

Xjj denotes a lattice site (i,j). The array iocc keeps track of the occupation 
status of the sites of the square lattice. If the site x^ is occupied then 
iocc(i,j)—l, otherwise it is equal to zero. Eqn.(2) expresses the fact that 
diffusion to a site from a neighbouring site takes place only if the neighbouring 
site belongs to the RD network. One can easily check that the discretized 
differential equations (with = \i a and ph = fJ-h) have a steady state 

solution given by a=l and h=l for all the lattice sites. Random fluctuations 
of magnitude less than 0.1 are created in the steady state with the help of the 
random number generator. This fixes the values of a and h at all the cluster 
sites at time t=0. The values of a and h at time t+1 are determined at a site 
Xij belonging to a cluster with the help of the discrete equations for a and 
h. This process is repeated till the steady state is reached, i.e., the values 
of a and h at all the cluster sites do not change within a specified accuracy. 
We define an 'activated' site as one at which the activator concentration has 
a value greater than 1, which is the homogeneous steady-state value. Figs. 
l(a)-l(c) show the patterns of activated sites on a lattice of size 50x50 and 
for site occuptation probability p = 1.0, p = 0.9, and p = 0.59. Fig. 1(a) 
shows a fully-connected stripe pattern with long-range connectivity. For p 
greater than a critical value p c , a cluster of n.n. activated sites still spans 
the lattice. Fig. 1(b) shows a pattern in which long-range connectivity still 
exists. Fig. 1(c) shows a pattern without long-range connectivity (p < p c ). 
The spanning of the largest cluster is checked as in the percolation problem 
and the value of p c is measured by the method of binary chopping fl7fl . Table 
I shows the average p c values determined for lattices of various sizes. For each 
lattice size, p c is calculated by taking the average of 100 values. Because of 
the finite size of the lattice, the value of p c has a small spread. As in the 
percolation problem, one expects that p c has a sharp, unique value as the 
lattice size L— > ex. The critical value p c is greater than the site percolation 
threshold p^ tte =0.59 for the square lattice. 

We now describe some cluster-related properties. Fig. 2 shows the span- 
ning cluster for RD on a 50x50 lattice and at p=p c (Table I). Let M be the 
size of the spanning cluster, size denotes the number of sites in the cluster. 
Size of the largest cluster is determined for lattice size LxL, L ranging from 
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20 to 70 in steps 10. To calculate M for each lattice size, the average of 200 
values is taken. Table II shows the data for average M for various lattice 
sizes. The fractal dimension df of the largest cluster is defined as 

M ~ L d f (3) 

Fig. 3 shows a plot of log(M) vs. log(L). From the slope, the value of df 
is obtained as <iy=1.86± 0.01 . The error quoted is the least-square fitting 
error. The value of df is close to that of the spanning cluster (dj=1.89) of 
ordinary percolation. At any value of p<p c , the non-equilibrium steady-state 
pattern consists of clusters of various sizes. Let n s be the number of clusters 
of size s per site of the lattice. The average cluster size Sav is then defined 
as 

S A v~± s s 2 n s (4) 

The prime indicates that the largest cluster is not included in the sum. As 
p — >p c , Sav exhibits a power-law scaling of the form 

Sav~\p~ Pel" 7 (5) 

The critical exponent 7 can be determined by various numerical methods. 
We use the method of finite-size scaling to obtain an estimate of 7 in our 
model. In this method, the average cluster size as a function of p and L is 
given by 

S AV (p,L) ~ Li f((p-p c )L») (6) 

The exponent v is the correlation length exponent. As p— >p c , the correlation 
length £ diverges as 

f ~ \P-Pc\~" (7) 

The function / in Eqn.(6) is a suitable scaling function and at p=p c , /(0)=constant. 
The average cluster size is determined for lattice size L ranging from 20 to 
70 in steps of 10. For each lattice size, the average of 200 values is taken. 
Table III shows the data obtained for p=p c . Fig. 4 shows a plot of \og{S A v) 
versus log(L) for p=p c - From the slope of the straight line, the exponent - 
is estimated as ^ = 1.72 ± 0.01. Again, the error quoted is the least-square 
fitting error. The value of - = 1.79. The average cluster size Sav is further 
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determined for (p — p c ) in the range 0.03 to 0.33. Fig. 5 shows a plot of 



70) are shown using the symbols: solid circle (L=50), open circle (L=60) 
and solid square (L=70). According to Eqn. (6), all the data should fall on 
a single curve with the functional form f(z). The value of the correlation 
length exponent v for which the best collapse of data is obtained is v = 1.66. 
The value of v is v = 1.75 in the case of ordinary percolation. 

III. Conclusion 

We have studied a RD model on the square lattice. The non-equilibrium 
steady state of the model consists of stripe-like patterns with long-range 
connectivity. A stripe is a region of activator concentration greater than the 
homogeneous, steady state value of 1. The sites of the underlying square lat- 
tice are occupied with probability p. As p falls below a critical value p c , the 
long-range connectivity of the steady state pattern disappears. As p — > p c , 
the average cluster size diverges. The spanning cluster at p = p c is a fractal 
object. The values of the exponent ^ and the fractal dimension df appears 
to indicate that the percolation model belongs to the universality class of 
ordinary percolation. The value 1.66 of the correlation length exponent v is 
different from the value 1.75 in the case of ordinary percolation. However, the 
difference may be due to finite-size effects and also due to the fact that the 
number of trials (~ 200) over which an average is taken is not large. This is 
because each trial takes a considerable amount of time. First, a steady state 
has to be generated which is reached only after several thousands of time 
steps. It is for this very reason that the lattice size could not be made large. 
The surface reaction models, we have discussed before, belong to the univer- 
sality class of directed percolation whereas the GM model of RD appears to 
belong to the universality class of ordinary percolation. In the ordinary per- 
colation problem, one looks at the connectivity of sites which are randomly 
occupied with probability p. The occupation status of a particular site is 
independent of the occupation status of neighbouring sites. In our study, we 
focus on the connectivity of a RD steady state pattern formed on a disor- 
dered lattice. The sites of the pattern are not randomly occupied, as in the 
case of ordinary percolation, but there is some degree of correlation in the 
occupancy of sites. The distribution of sites in the steady-state depends on 




versus 




Data for three lattice sizes (L=50, 60 and 
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the parameters of the RD model, the diffusion coefficients D a , Dh, the sat- 
uration parameter k a etc. The cluster-related properties of the steady state 
pattern, say, the average cluster size, have a power-law scaling behaviour 
when occupation probability p of the underlying lattice is close to a critical 
value p c . As already mentioned in the Introduction, there are many physi- 
cal, chemical and biological systems which exhibit Turing patterns. Some of 
these patterns consist of stripes which may sometimes form in a disordered 
environment. The present study is of relevance to such systems. In the per- 
colation picture, one studies a transition from localised to extended states 
as a function of disorder. In the GM model, instead of disorder, one can 
vary the parameter k a to bring about the transition. For k a =0, the steady 
state pattern consists of isolated Turing islands. For k a ^ 0, the steady 
state consists of connected stripes. An interesting question to ask is whether 
the transition from localised to extended states occurs as soon as k a is made 



different from zero. Preliminary studies [IS shows that the transition occurs 
at a finite value of k a . The connection of this transition with the one in the 
presence of disorder will be explored in future. The effect of changing the 
parameters of the RD model and also changing the form of the nonlinearity 
in the RD model (provided it still gives rise to connected patterns in the 
steady state) will also be investigated. 
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Table I. The average value of p c for lattices of various sizes. 



Latice size 


Average p c 


20 x 20 


0.81089 


30 x 30 


0.80735 


40 x 40 


0.80992 


50 x 50 


0.80714 


60 x 60 


0.80873 


70 x 70 


0.80843 



Table II. The average size M of the spanning cluster at p = p c for lattices of 

various sizes 



Latice size 


Average M 


20 x 20 


96.54 


30 x 30 


195.77 


40 x 40 


350.65 


50 x 50 


508.77 


60 x 60 


754.81 


70 x 70 


974.22 



Table III. The average cluster size Sav at V — Pc for lattices of various sizes 



Latice size 


Average Sav 


20 x 20 


4.01 


30 x 30 


8.22 


40 x 40 


14.11 


50 x 50 


19.65 


60 x 60 


27.82 


70 x 70 


33.96 
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Figure Captions 

Fig.l Non-equlibrium steady state patterns in the GM model on a disor- 
dered square lattice of size 50 x 50. The probability that a site of the 
original square lattice exists is p. The patterns are shown for (a) p—1, 
(b) p=0.9, (c) p=0.59 . 

Fig. 2 The spanning cluster on a 50 x 50 lattice for p = p c (see Table I.) . 

Fig. 3 log(M) versus log(L) for p = p c (Eqn.(3)). The slope of the straight 
line gives the fractal dimension rf/=1.86±0.01 . 

Fig. 4 \og(SAv) versus log(L) for p = p c (Eqn.(6)). The slope of the straight 
line gives the exponent ^=1.72±0.01 . 

Fig. 5 The collapse of data for three different lattice sizes: solid circles, open 
circles and solid squares for L=50, 60 and 70 respectively. 
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